Consistency Analysis of Finite Difference 
Approximations to PDE Systems 



Vladimir P. Gcrdt 

Laboratory of Information Technologies, 

Joint Institute for Nuclear Research, 141980 Dubna, Russia 

gerdt© . j inr . ru 



(~~) • Abstract. In the given paper we consider finite difference approxima- 

tes) ' tions to systems of polynomially-nonlinear partial differential equations 

whose coefficients are rational functions over rationals in the indepen- 
O ■ dent variables. The notion of strong consistency which we introduced 

^^ ' earlier for linear systems is extended to nonlinear ones. For orthogonal 

and uniform grids we describe an algorithmic procedure for verification 
of strong consistency based on computation of difference standard bases. 
The concepts and algorithmic methods of the present paper are illus- 
trated by two finite difference approximations to the two-dimensional 
Qh ' Navier-Stokes equations. One of these approximations is strongly consis- 

tent and another is not. 
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*0 ■ 1 Introduction 

> 

v^ ' Along with the methods of finite volumes and finite elements, the finite difference 

f»o , method ^24 is widely used for numerical solving of partial differential equations 

(PDE). This method is based upon the application of a local Taylor expansion 
r — | , to replace a differential equation by the difference one [26128] defined on the cho- 

f^ ' sen computational grid. The last equation forms finite difference approximation 

(FDA) to the given PDE, and together with discrete approximation of initial 
or/and boundary condition constitutes a finite difference scheme (FDS). 

In theory, the most essential feature required of discretization is convergence 
of a solution of FDS to a solution of PDE as the grid spacings go to zero. However, 
except a very limited class of problems, convergence cannot be directly analyzed. 
Instead, it has been universally adopted that convergence is provided if FDA 
is consistent and stable. This adoption is due to the brilliant Lax-Richtmyer 
equivalence theorem [26 28 proved first for linear scalar PDE equations and 
then extended to some nonlinear scalar equations |23| . The theorem states that 
a consistent FDA to a PDE with the well-posed initial value (Cauchy) problem 
converges if and only if it is stable. Consistency implies reduction of FDA to the 
original PDE when the grid spacings go to zero. It is obvious that consistency is 
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necessary for convergence. As to stability, it provides boundedness of the error 
in the solution under small perturbation in the numerical data. 

Thus, the consistency check and verification of stability are principal steps in 
qualitative analysis of FDA to PDE. Modern computer algebra methods, algo- 
rithms and software may provide a powerful tool for generating FDA [11] and for 
performing its consistency and stability analysis. Some recent computer algebra 
application to study stability and to generate FDA to linear PDE systems with 
constant coefficients are discussed in (50] • In papers [10113] some computer alge- 
bra and algorithmic issues related to the consistency analysis were considered. In 
particular, for orthogonal and uniform solution grids the notion of s-consistency 
(strong-consistency) was introduced in [T3] for FDA to a linear PDE system 
that strengthens the conventional notion of consistency and admits algorithmic 
verification. In doing so, an s-consistent discretization not only approximates the 
differential equations in a given linear system but also preserves at the discrete 
level algebraic properties of the system. It follows that if the system has local 
conservation laws in the form of algebraic consequences of its equations, then 
the s-consistent discrete system will also have such conservation laws (cf. |5I29] V 

In this paper we generalize the concept of s-consistency to polynomially- 
nonlinear PDE systems and extend the algorithmic ideas of paper 13 to check 
s-consistency for such systems on orthogonal and uniform solution grids. In the 
linear case algorithmic verification of s-consistency is based on completion of the 
initial differential system to involution and on construction of a Grobner basis for 
the linear difference ideal generated by FDA. It is important to emphasize that 
involutivity of the linear differential system under consideration not only makes 
possible an algorithmic verification of s-consistency but is also necessary (cf. |25j ) 
to well-posedness of Cauchy problem for the system what, if one believes in the 
extension of Lax-Richtmyer equivalence theorem to PDE systems, can provide 
convergence for s-consistent and stable FDA. 

However, a differential system may not admit involutive form. Generally, 
one can decompose such a system into a finitely many involutive subsystems by 
applying the Thomas decomposition method |27) . The decomposition is done 
fully algorithmically [TJ with the use of constructive ideas by Janet [TB] further 
developed and generalized in |7I9| . Another obstacle for nonlinear FDA is that 
the relevant nonlinear difference Grobner basis |19| may be infinite. Since it is 
commonly supposed that Grobner basis is a finite object, its infinite difference 
analogue is called standard basis as well as in differential algebra (cf. [21131] ). 

This paper is organized as follows. Section 2 contains a short description of 
differential and difference systems of equations which are studied in the paper. 
The properties of differential Thomas decomposition that are used for the s- 
consistency check are considered in Section 3. In Section 4 we define difference 
standard bases and present an algorithm for their construction. The definition of 
s-consistency of FDA for uniform and orthogonal grids, which is a generalization 
of that in [T5] to nonlinear differential systems, is given in Section 5. Here we 
also formulate and prove the main theorem on the algorithmic characterization 
of s-consistency and propose an algorithmic procedure for its verification. The 
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concepts and methods of the paper are illustrated in Section 6 by two FDA 
derived in [TU] for the two-dimensional Navier-Stokes equations. Some concluding 
remarks are given in Section 7. 

2 Preliminaries 

In the given paper we consider PDE systems of the form 

.fi = --- = f P = 0, F:={f u ...J p }cK. (1) 

Here /, (i — 1, . . . ,p) are elements in the differential polynomial ring TZ := 
IClu 1 , . . . , u m ], that is, polynomials in the dependent variables u := {u 1 , . . . , u m } 
{differential indeterminates) and their partial derivatives which are the operator 
power products of the derivation operators {Si, . . . ,5 n } (Sj = d x ). We shall 
assume that coefficients of the polynomials are rational functions in the inde- 
pendent variables x := {x\, . . . , x n } whose coefficients are rational numbers, i.e. 
£:=Q(x). 

To approximate the differential system (HJ by a difference system we shall 
use an orthogonal and uniform computational grid (mesh) as the set of points 
(kihi, . . . ,k n h n ) in W 1 . Here h := (hi, ■ ■ ■ , h n ) (hi > 0) is the set of mesh 
steps (grid spacings) and the integer- valued vector (ki, . . . , k n ) G Z™ numerates 
the grid points. If the actual solution to the problem (|TJ) is the vector-function 
u(x), then its approximation in the grid nodes will be given by the grid (vector) 
function u^,...,^ = u(kihi,. . .,k n h n ). 

We shall assume that coefficients of the differential polynomials in F do not 
vanish in the grid points. The coefficients on the grid as rational functions in 
{fci/ii, . . . , k n h n } are elements of the difference field |1!J with mutually commut- 
ing differences {o~i, . . . , a n } acting on a function <f>(x) as the right-shift operators 

Oi O (j)(xi,. . .,£„) = 4>(Xl,. . . ,Xi + hi, . . . ,x n ) (hi>0). (2) 

The monoid (free commutative semigroup) generated by a will be denoted 
by 0, i.e. 

0:={o-?o...oo-i" |»i,...,t n eN>o}, (V6eO) [^ol = l], 



the field of rational functions in {kih\, . . . , k n h n } by /C and the ring of difference 
polynomials over /C by 1Z. The elements in 1Z are polynomials in the dependent 
variables (difference indeterminates) u a (a — 1, . . . , m) defined on the grid and 
in their shifted values a 1 ^ o- • -ocr*™ ou" (ij g N>o). The coefficients of polynomials 
are taken from K,. 

The standard technique to obtain FDA to ((T|) is replacement of the deriva- 
tives occurring in (JT|) by finite differences and application of appropriate power 
product of the right-shift operators @ to remove negative shifts in indices which 
may come out of expressions like 

W _ M 

„ (i) u ki,...,kj+l,...,k n u k ls ...,k j -l,...,k n /n / ; ,2\ 

dj o u w = — h O(hj) . 
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In [TT] we suggested another approach to generation of FDA which is based 
on the finite volume method and on difference elimination. As it was shown 
for the classical Falkowich-Karman equation in gas dynamics, this method may 
derive a FDA which reveals better numerical behavior then those obtained by 
the standard technique. In the sequel we shall consider FDA to the PDE system 
(fT]) as a finite set of difference polynomials 

/i = --- = /, = 0, F:={h,...J q }cH, (3) 

where q need not be equal to p. 

We shall say that a differential (resp. difference) polynomial / € H (resp. 
/ G TV) is differential- algebraic (resp. difference- algebraic) consequence of (fTJ) 
(resp. ([3])) if / (resp. /) vanishes on the common solutions of (P) (resp. ©). 



3 Differential Thomas Decomposition 

Definition 1. Let S = and S^ be finite sets of differential polynomials such 
that S = ^ and contains equations (Vs e S = ) [s — 0] whereas S^ contains 
inequations (Vs £ .fv ) [s / 0]. Then the pair [S = ,S^) of sets S = and <fv is 
called differential system. 

Let &ol(S = / S^) denote the solution set of system (S* = ,5^), i.e. the set 
of common solutions of differential equations {s = | s 6 5 = } that do not 
annihilate elements s 6 S^ . 

Theorem 1. |27j Any differential system [S = , S*^) is decomposable into a finite 
set of involutive subsystems yS^~,Sf) with disjoint set of solutions 

(s=/s*) =► |J (s=/sf) , Goi(sr/s*) = |+| eoi (sr/sf) . (4) 



The structure of involutive subsystems in decomposition of a given system 
depends on the choice of ranking defined as follows. Consider the monoid of 
derivation operators A :— { S^ o ■ • • o S l ™ \ii,.. .,i n € N>o }• 

Definition 2. A total (linear) ordering y on the set of partial derivatives {So 
u a | 8 € A, a = 1, . . . ,p} is ranking if for all i, a, /?, 6, 5 

5i o Su a y S o u a , 5 o u a >- 5 o u 13 •<=>• 5i o 6 o u a y Si o 5 o u@ . 

If (3 7 ) [Sou 1 y Souf] => (Va,/3) [Sou a y Sou 13 }, then y is orderly. // 
u a y u@ =>■ (V S,S) [5 o u a y 5 o u@ ], then y is elimination. 

The Thomas decomposition into Janet involutive |16j subsystems is done 
fully algorithmically and have been implemented as a Maple package Q]. Given 
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decomposition Q , one can algorithmically verify whether a differential equation 
/ = (/ G TZ) is a differential-algebraic consequence of the system (S = , S^) 

(Va £ Got (S=/S*) [/(o) = 0]^^(Vi)[d P rem i7 (/ ) ST)=0]. (5) 

Here dprenij-(/, P) denotes differential Janet pseudo-reminder of f modulo P. 
The underlying Janet pseudo-division algorithm is described in [T] and imple- 
mented in the package. 

Remark 1. For the case S^ = condition ([5]) verifies / G [<S = ] C TZ, where 
[-F] denotes the radical of differential ideal generated by the set F. Thereby, the 
Thomas decomposition of (F, 0) provides a characteristic decomposition of |F] 
(see [1114) for more details). 

Example 1. We illustrate the Thomas decomposition by the the example taken 
from [5]. Consider differential system 

({(uy + v)u x + 4v u y — 2v , (u y + 2u)u a; + 5v u y — 2v }, {}) 

with two quadratically-nonlinear first-order PDE with two dependent and two 
independent variables. Its Thomas decomposition for the ranking satisfying u x y 
u y y v x y v v y u y V is given by 

(u y + v)u x + 4v u y — 2v 2 \ , V , 

For a differential system with linear PDEs and the empty set of inequations 
the decomposition algorithm performs completion of the system to involution 
and returns the Janet basis form [316] of the input system. 

4 Difference Standard Bases 

For the shifted dependent variables ranking is defined in perfect analogy to Def- 
inition [5] of ranking for partial derivatives. 

Definition 3. \F3\ A total ordering -< on {8 o u a | 9 G 6, 1 < a < m} is 
ranking if for all <jj , 9, 9\ , 9i , a, /3 

(i) a, o 9 o u a y 9 o u a , (ii) 6 X o u a y 9 2 o ?/ ^=^ o 6>i o u a y 9 o 6» 2 o i/ . 



Definition 4. ^4 totaZ ordering y on the set Ai of difference monomials 
M := { (6»i o u 1 ) 11 ■ ■ ■ (9 m o u m y™ | 9j G O, i 3 G N> , 1 < j < m } 
is admissible z/ «i extends a ranking and satisfies 
(VteM\{l}) [tyl] A (V0e6>) (Vt,u,weM) [w >- tu •$=>• M9o« >- t-9o W }. 
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Remark 2. Similar to that in Definition [2] one can define orderly and elimination 
difference rankings. As an example of admissible monomial ordering we indicate 
the lexicographical monomial ordering compatible with a ranking. 

Given an admissible ordering >-, every difference polynomial / has the leading 
monomial lm(/) 6 M. with the leading coefficient lc(/). In what follows every 
difference monomial is to be normalized (i. e. monic) by division of the monomial 
by its leading coefficient. This provides (V/ £ TZ) [lc(/) = f ]. 

Now we consider the notions of difference ideal [TH] and its standard basis. 
The last notion is given here in analogy to that in differential algebra [21] . 

Definition 5. [19] A set I C TZ is difference polynomial ideal or cr-ideal if 

(Va,6el) (Vceft), (V6>£<9) [o + kl, a-cel, 6oael]. 

If F C TZ, then the smallest a-ideal containing F is said to be generated by F 
and denoted by [F] . 

If for v, w £ M. the equality w = t ■ 9 o v holds with 9 £ and t £ M we 
shall say that v divides w and write v | w. It is easy to see that this divisibility 
relation yields a partial order. 

Definition 6. Given a a-ideal I and an admissible monomial ordering >-, a 
subset G <zX is its (difference) standard basis if [G] = X and 

(Vfel)(3~geG) [lm(.9)|lm(/)]. (6) 

If the standard basis is finite it is called Grobner basis. 

Definition 7. A polynomial p £ TZ is said to be head reducible modulo q £ 1Z 
to r if f = p — m ■ 9 o q and m £ A4, 9 £ O are such that lm(j>) = m ■ 9 o hn(g). 
In this case transformation from p to f is elementary reduction and denoted by 
p — y f. Given a set F C TZ, p is head reducible modulo F (denotation: p — >) if 

? F 

there is f £ F such that p is head reducible modulo f . A polynomial p is head 
reducible to f modulo F if there is a chain of elementary reductions 

p — >pi — >p 2 —>•■•■—> f. (7) 

F F F F 

Similarly, one can define tail reduction. Iff in (p|) and each of its monomials are 
not reducible modulo F , then we shall say that f is in the normal form modulo 
F and write f = NF(p, F). A polynomial set F with more then one element is 
interreduced if 

(V/eF)[/=NF(/,F\ {/})]. (8) 
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Admissibility of >-, as in commutative algebra, provides termination of chain ([7]) 
for any p and F. In doing so, NF(p, F) can be computed by the difference version 
of a multivariate polynomial division algorithm [214] . If G is a standard basis of 
[G] , then from Definitions [6] and [Jj it follows 

/G[G]^NF(/,G) = 0. 

Thus, if an ideal has a finite standard (Grobner) basis, then its construction 
solves the ideal membership problem as well as in commutative [214] and differ- 
ential [21131] algebra. The algorithmic characterization of standard bases, and 
their construction in difference polynomial rings is done in terms of difference 
S'-polynomials. 

Definition 8. Given an admissible ordering, and monic difference polynomials 
p and q, the polynomial S(p, q) := m\ ■ 9\ op — mi ■ #2 ° q is called S'-polynomial 
associated to p and q (for p = q we shall say that S-polynomial is associated 
with p) if mi ■ 0\ o lm(j>) = m^ • #2 ° lm(g) with co-prime wi\ ■ B\ and m-i • Qi- 

Theorem 2. Given an ideal I C R and an admissible ordering )~, a set of 
polynomials G CI I is a standard basis of I if and only j/NF(S(p, q), G) = for 
all S -polynomials, associated with polynomials in G. 

Proof. It follows from Definitions [BJ [7] and |5] in line with the standard proof of 
the analogous theorem for Grobner bases in commutative algebra |2I4) and with 
the proof of similar theorem for standard bases in differential algebra [2T] . □ 

Let I — [F] be a cr-ideal generated by a finite set F C 1Z of difference poly- 
nomials. Then for a fixed admissible monomial ordering the below algorithm 
StandardBasis, if it terminates, returns a standard basis G of X. The subal- 
gorithm Interreduce invoked in line 11 performs mutual interreduction of the 
elements in H and returns a set satisfying ((8]). 

Algorithm StandardBasis is a difference analogue of the simplest version 
of Buchberger's algorithm (cf. |2l4l21j h Its correctness is provided by Theorem 
[2] The algorithm always terminates when the input polynomials are linear. If 
this is not the case, the algorithm may not terminate. This means that the do 
while- loop (lines 2-10) may be infinite as in the differential case [21131] , One 
can improve the algorithm by taking into account Buchberger's criteria to avoid 
some useless zero reductions of line 5. The difference criteria are similar to the 
differential ones [2T] . 

Example 2. Consider a simple example of the principal ideal generated by poly- 
nomial gi :— u(x) ■ u(x + 2) — x ■ u(x + 1) in the ordinary difference ring with 
the only shift operator (difference) aou(x) — u(x + 1), the independent variable 
(indeterminate) u and the dependent variable x. Let us fix monomial ordering 
as the pure lexicographic one with u(x) -< u(x + 1) -< • • • . Obviously, it is ad- 
missible. Then a nontrivial (i.e. having nonzero normal form) S'-polynomial s\ 
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associated with g\ and its normal form 32 modulo {31} are given by 

Si := u(x + 4) • .gi - u(x) ■ a 2 o g ± , 

.t + 2 

9 2 := NF(si, {<h}) = m(.x + 1) • u(x + 4) - ■ u(x) • u(a; + 3) . 

x 

The second nontrivial S'-polynomial S2 associated with g\ , g 2 and its normal form 
g 3 modulo {51,32} read 

s 2 := u(x + 4) • <j o g x — u(x + 3) • g 2 , 

g 3 := NF(s 2 , {31,32}) = u(x) ■ u(x + 3) 2 - x ■ (x + 1) • u(x + 3) . 

One more nontrivial ^-polynomial S3 associated with 32 , 33 and its normal form 

34 modulo {31,32,33} are 

S3 := cr o -33 - u(x + 4) • 32 , 

34 := NF(s 3 , {31, 32, 33}) = u(x) ■ u(x + 3) • u(x + 4) - x • (a; + 1) • u(x + 4) . 

The last nontrivial S'-polynomial S4 associated with 33,34 and its normal form 

35 modulo {31,32,33,34} are 

s 4 := u(x + 5) • 33 - cr o 34, , 

x + 3 

35 := NF(s 4 , {31,32,33,34}) = u(x + 5) — r _ 7T' u ( a; ) ' u ( x + 4 ) ■ 

a; • (x + 1) 

Now all 5-polynomials associated with elements in G := {31,32,33,34,35} are 
reduced to zero modulo G, and G is an interreduced standard basis of [31]. 

Algorithm: StandardBasis (F , >-) 



Input: F € TZ\ {0}, a finite set of nonzero polynomials; 

>-, a monomial ordering 
Output: G, an interreduced standard basis of [F] 



G:=F 
do 

i7:=G 

for all S'-polynomials s associated with elements in H do 
3 :=NF(s,#) 
if 3 ^ then 

G:-GU{3} 
fi 
od 
od while G^il 
G :=Interreduce (G) 
return G 
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5 Consistency of Finite Difference Approximations 

For simplicity, throughout this section we shall consider orthogonal and uniform 
grids with equisized mesh steps h± = ■ ■ ■ = h n = h. 

Definition 9. [T3] We shall say that a difference equation /(u) = implies the 
differential equation /(u) = and write f > f when the Taylor expansion about 
a grid point yields 

/(u) ► f(u)h k + 0(h k+1 ), k e Z >0 . 

h— >0 

Definition 10. (13) Given a PDE system (QJ) and its difference approximation 
(0), we shall say that |[3J) is weakly consistent or w-consistent with (QJ) if 

(V/eF) (3feF) [/>/]. 



In paper [T3] we showed that already for linear PDE systems such definition of 
consistency, which has been universally accepted in the literature, is not sat- 
isfactory in view of inheritance of properties of differential systems by their 
discretization. Instead, we introduced another concept of consistency for linear 
FDA which is extended to nonlinear systems of PDE as follows. 

Definition 11. [TI5] A perfect difference ideal generated by set F 6 TZ and de- 
noted by \F\ is the smallest difference ideal containing F and such that for any 
f e U, 0i,..., r e and ki,...,k r € N> 

(6 X o /> ■■■(9 r o /> e {F} =► / e [F] . 

It is clear that [F] C [F]. In difference algebra perfect ideals play the 
same role as radical ideals in commutative [3] and differential algebra [14 , for 
example, in Nullstellensatz [30] . By this reason we shall consider the perfect 
ideal \F\ generated by the difference polynomials in FDA <j3j> as the set of its 
difference-algebraic consequences. Respectively, the set of differential- algebraic 
consequences of a PDE system is the radical differential ideal generated by the 
set F in ([1]) (see Remark [TJ. 

Definition 12. An FDA £3$) to a PDE system UJ) is strongly consistent or 
s-consistent if 

(V/e[F]) (3/ e[FJ) [/>/]. (9) 



The algorithm ConsistencyCheck presented below verifies s-consistency of 
FDA to PDE systems. Its correction is provided by property ([5]) of the differential 
Thomas decomposition applied in lines 13-16 of the algorithm and by Theorem 
[3] This theorem generalizes to nonlinear systems the theorem formulated and 
proved in 113: for linear systems. 



10 Vladimir P. Gerdt 

Theorem 3. A difference approximation (0) to a differential system (QJ) is s- 
consistent if and only if a reduced standard basis G C 1Z of the difference ideal 
[F] satisfies 

(Vg€G)(3gem)[g>g]. (10) 



Proof. Let >- be an admissible monomial ordering and G be the correspond- 
ing interreduced standard basis. To prove that (|10l) implies ^ consider first a 
nonzero polynomial / G [F] and show that / D> / G [FJ . Polynomial / as well as 
any ^-polynomial associated with elements in G, because of the property ([B]) of 
G, admits representation with respect to G and >- as a finite sum 

/= Yl '^2 a 9^- (J ' 1 °9, ag tli €%, Im(a §)/i ■ a* o g) < lm(/) . (11) 

gGGiCG M 

Here we use the multiindex notation 

H := 0*1, • • • , /V) e Z| , <r" := of o • • • o ct£" . 

Formula dill) is a difference analogue of the standard representation in com- 
mutative algebra [2|. Consider the Taylor expansion (in grid spacing h) of the 
right-hand side of (TlTT) about a grid point, nonsingular for the coefficients occur- 
ring in the sum. In doing so, the shift operators aj (j = 1, . . . , n) are expanded 
in the Taylor series 

fc>0 

along with the shifted coefficients as rational functions in the independent vari- 
ables. 

The representation (|llj) guarantees that in the leading order in h the leading 
differential monomials |21j which occur in the sum and come from different 
elements of the Grobner basis cannot be cancelled out. Thereby, due to the 
condition (|10p , the Taylor expansion of / implies a finite sum of the form 

/ : =EE^-^°fl, b g , v &n, 

g£Gi (U 

where d := {g <=Tl \ 3g E G x such that g > g}. Therefore, / > / e [F] C {Fj. 
Let now ft G [F] \ [F] and 9\ , . . . , 9 r £ and k\ , . . . , k r e N> be such that 

9:=(^°^'''^op) fer e[f]. (13) 

As we have shown, q\>q G [F], and it follows from (|P^|) that g = p fclH hfcr where 
pt>p. Hence, p G \F\. The perfect ideal [F] can be constructed [TO] from [F] 
by the procedure called shuffling and based on enlargement of the generator set 
F with all polynomials p satisfying (1131) and on repetition of such enlargement. 
It is clear that each such enlargement of the intermediate ideals yields in the 
continuous limit a subset of [F]. 

Conversely, conditions (|10l) trivially follow from © and from G C [FJ. □ 
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Algorithm: ConsistencyCheck (F, F) 





Input: F C TZ \ {0}, F £ 1Z \ {0}, finite sets of nonzero polynomials 


Output: true if F is s-consistent FDA to F, and false otherwise 


1 


choose differential ranking >~i and difference ordering ^2 


2 


T :—DifFerentialThomasDecomposition(F, Xi) 


3 


V :={P\(P,Q)eTl 


4 


G := StandardBasis (F, ^2) (* may not terminate *) 


5 


C := true 


6 


while G 7^ and C — true do 


7 


choose q £ G 


8 


G:=G\{g}; V := V 


9 


compute g such that g\> g 


10 


while V 7^ and C = true do 


11 


choose S E P 


12 


P :=P\{S} 


13 


d:=dpremj(g,S) 


14 


if d ^ then 


15 


C := false 


16 


fi 


17 


od 


18 


od 


19 


return C 



ft should be noted that condition (O does not exploit the equality of cardinal- 
ities for sets of differential and difference equations as is assumed in Definition 
[TU1 The equality of cardinalities is also not used in the proof of Theorem [31 
Therefore, both Definition [T21 and Theorem [3] are relevant to the case when the 
FDA has the number of equations different from that in the PDE system. 

In the nonlinear case when algorithm StandardBasis may not terminate, 
it is useful to compute the continuous limit g > g for the difference polynomials 
g obtained in line 5 of algorithm StandardBasis and to verify the condition 
dpremj(g, S) — as it is done in line 14 of algorithm ConsistencyCheck. 
This way one can stop computation when inconsistency of the intermedite data 
in algorithm StandardBasis is detected. An example of such situation is con- 
sidered in the next section. 



6 Example: Navier-Stokes Equations 

To illustrate the concept of s-consistency and the algorithmic procedure of its 
verification, we consider two FDA generated in [10) for the two-dimensional 
Navier-Stokes equations by the method of paper [TT]. These equations describe 
unsteady motion of incompressible viscous liquid of constant viscosity. The Janet 
involutive form of the Navier-Stokes equations for the orderly ranking compatible 
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with S x >~ S y >~ St and u >~ v >~ p is given by (see jlOj ) 



F := 



h 


— U X + Vy = , 




h 


= U t + UU X + VUy +P X - -^(U xx + Uyy) 


= 


h 


= V t + UV X + VVy + Py - -^{V XX + Vyy) = 


= 0, 


fa 


= U\ + 2V x Uy + l'l + P XX + Pyy = . 





(14) 



Here fa is the continuity equation, fa and fa are the proper Navier-Stokes equa- 
tions 22], fa the pressure Poisson equation [T3], (u, v) is the velocity field, and 
p is the pressure. The density is included in the Reynolds number Re. 

The differential Thomas decomposition algorithm [T] for the input fi, fa, fa 
outputs system fTJ| in its Janet autoreduced form 

U x + Vy = , 

r i HiKi/ ~ v *y ~ uv v) -vu v -ut-p x =0, 

RSV^z + %j/) _ uv x ~ VV V ~Vt-Py = 0, 

2V X Uy + P XX + Pyy + 2Vy = . 

The following FDA to system (fT4|) was obtained in [10] for the orthogonal 
and uniform grid with the spatial spacing h and temporal spacing r: 



/l — 2h 1 2h - U ' 



fa 



ti * "3 fc I " j + lfc~" j-lfc I ""jfc + 1- 



2h ' 2/i 



Pj + l fe Pj-l fe 
2to 



J_ f u j + 2k~ 2u "k+ u "-2k 1 ""fc + 2~ 2 ""fc+""fc-2 N _ r. 

Re \ IP ' iH 2 J ~ U • 



? ._ v 7i 1 - v ik , "^"+ifc-"«"-ifc , -" 2 J " fc +i-^ 2 " fc -i , 

« '~ t T 2h " r 2/i "f" 

"T" 2h Re ^ UP ' UP J ~ U • 

f u2 "+2k- 2u2 "k+ u2 ]-2k 1 r, ut '"+lfe + l- t "'"+lfc-l- ut '"-lfc + l+ Mt '"-lfc-l 

J4 •- jui h ^ — ■ 



4/i 2 : - 4h 2 

,-2i> 2 



.7 fc ' u jfe — 2 



4/i 2 



I ( p]+2k~ 2 P7k+P?-2k 1 P"fc + 2- 2 P"fc+P"fc-2 A 

T ^ 4/i 2 + 4/i 2 y 



= 0. 



This FDA is w-consistent what can be easily verified by the Taylor expansion 
of the difference polynomials in F :— {fa, fa, fa, fa} in the powers of h, r about 
a grid point. In doing so, in the continuous limit (r — > 0, /i — ► 0) the difference 
equations imply the involutive differential Navier-Stokes system (fT4"l) . Moreover, 
the algorithm StandardBasis applied to the set F\ :— {<j y o fa,o y o fa, fe, fi\ 
yields that F\ is a difference Grobner basis of ideal [F{\ for the lexicographic 
ordering compatible with the orderly ranking such that a t >- a x >- a y and 
p>- u>- v (see Remark (2). Thus, F is the s-consistent FDA to (IT4^I . 

The above given FDA has a 5 x 5 stencil owing to the approximation of 
the second-order partial derivatives used for equations /2,/3 and fa. From the 
numerical standpoint a 3 x 3 stencil looks like more attractive. By this reason 
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let us follow [TO] and consider another FDA to (fT4")l with a 3 x 3 stencil: 



2 "j + l k "j-lfc , "jfc + 1 "jfc-l _ r, 

C l • Oh 9h u • 



s ._ J jfc -" 3 -fc , " j+ifc-" j-ifc , uv j k+i-™ j k-i , 

c 2 ■— T -T- 2h "•" 2/i " r 

i Pj+ifc-Pj-ifc __ 1 ( u "+lk- 2u ™k+ u ™-lk , M "fc + l- 2 ""fc+""fc-l A _ n 
"•" 2/i Re ^ /i 2 "T" /i 2 / ' 

= ._ ^"fc 1 -""^ , ™?+ifc- ti '" J "-ifc , ^ 2 ?fc+i-" 2 7fc-i I 

e 3 •— T T 2 /i I" 2/> " r 

i p"fc+i-p"fc-i _ i /' t, "+ifc~ 2t '"fc+ t '"-ifc i •" J "fc+i- 2- ""fc+ , "7fc-i A _ n 

"*" 2/i Re ^ /i 2 ~r /i 2 / ' 

s '» 2 "+lfc- 2M2 "fc+" 2 "-lfc , ™-'j + ifc + i-™? + ifc-i-™"-ifc + i+™"-ifc-i i 

e 4 ■- P I" z iTT^ r 

i " jfc+i- 2 " jfc+" Jfc -i ^ p"+i fc~ 2 p"fc+p"-i ), , p"fc+i~ 2 p"fc+y"fc-i \ _ a 
H T? v \ 7? ! S 2 J~ u 

f := {ei, e2, e3, £4} is w-consistent with (|14j) . However, application of algo- 
rithm StandardBasis shows that, as opposed to F\ , F[ :— {<r y oei , a y oe.2, £3, £.4} 
it not a Grobner basis. For the S'-polynomial §1,2 associated with o y o §1 and 
Oyoe-i we have q :— NF(si i2 ,-Fi) 7^ 0. Furthermore, qOq := u xx + Vy y +p xx +p yy . 
The equation q — is not a consequence of the Navier-Stokes system. 

One way to check it is to compute d :=dpremj(q, F\) with F\ given by (fT5j) . 
Just this computation is done in line 13 of algorithm ConsistencyCheck: 

d = rV ( u yy + v yy ~ 2u y v * ~ 2v l) + Bi ( uv y u yy ~ vu y u yy ~ u t u yy - p x U w ) + 
2 {vu t u y — uu t v y + vu y p x - uv y p x — uWyU y + u t p x ) + u 2 + p x + w 2 ^ 2 , + u 2 u^ ■ 

Another way is to substitute into g the exact solution [IT] to (|14p 

u = — e~ 2t cos(x)sin(y), u = e~ 2 * sin(a;) cos(y), p= — e~ *(cos(2a;)4-cos(2j/))/4. 

and to see that it does not satisfy q — 0. Therefore, F' is s- inconsistent. 

7 Conclusion 

Our computer experiments |13) with linear systems based on the implementa- 
tion [TJI of Janet completion algorithm for the <r-ideals generated by linear dif- 
ference polynomials shown that unlike w-consistency it is fairly difficult to satisfy 
s-consistency by discretizing overdetermined PDE systems. This is hardly sur- 
prising since an s-consistent FDA preserves at the discrete level all consequences 
of the differential system. As we demonstrate in Section 6 of the present paper, 
completion of the Navier-Stokes equations to involution by adding the Poisson 
pressure equation, which has to be explicitly taken into account in the numerical 
solving [T3], makes the s-consistency of their FDA sensitive to discretization. 

To guarantee termination of the algorithmic versification of s-consistency one 
might use the fact that the difference polynomial ring we deal with in this paper 
is a Ritt ring and each its perfect ideal has a finite basis [TH]. However, unlike the 
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differential Ritt rings [14] , there are no algorithms known to compute such basis 
and, hence, a Grobner basis for [F]. Another obstacle in computer application 
to the consistency analysis of FDA is the lack of software for construction of 
nonlinear standard bases. Only very recently a start has been made with a new 
algorithmic insight inspired by the ideas of paper [TB] with intention to create 
such software packages written in Maple and Singula:^. 
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